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00 ■ Abstract 
0\ 

, Correlations among identical bosons, which are familiar from sta- 

Qh| tistical physics, play an increasingly important role in high energy mul- 

O ■ tiple particle production processes. They provide information about 

<D the region, where the particles are produced and, if Einstein's conden- 

sation can be reached, they can lead to spectacular new phenomena. 



h ■ 1 INTRODUCTION 



In this paper we will consider Bose-Einstein correlations in high-energy par- 
ticle production processes, i.e. the correlations among identical bosons in 
the final state, which follow from Bose-Einstein statistics. When hundreds 
of identical bosons are being produced in a single scattering act, as hap- 
pens e.g. in heavy ion collisions at high energy, such correlations can lead 
to spectacular phenomena. They are also, most probably, the best way of 
getting information about the space-time structure of the region, where the 
final state particles are produced. Let us begin with a very simple example. 
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Consider the elastic scattering of two alpha particles with initial momenta 
equal in magnitude, opposite and parallel to a horizontal axis, say the x-axis. 
Suppose that the detectors register the final state particles if and only if the 
scattering is at 90° and one of the final particles goes up and hits the upper 
detector (U), while the other goes down and hits the lower detector (L). 
There are two possibilities. Either particle 1, say the particle coming from 
the left, hits detector U and particle 2 hits detector L, or particle 1 hits 
detector L and particle 2 hits detector U. Let us denote the probability 
amplitudes for these two processes by A and B respectively. Since a rotation 
around the x-axis can convert these two processes into each other, \A\ = \B\. 
Since the alpha particles are identical bosons and an exchange of the two 
final state particles converts the two processes into each other, A = B. The 
detection probability is |A+_B| 2 = A\A\ 2 . If the particles were distinguishable, 
the probability would be \A\ 2 + \B\ 2 = 2\A\ 2 . Thus, the fact that the particles 
are indistinguishable increases the probability by a factor of two. 

Let us make a few comments about this simple result. 

• In the example the two amplitudes interfere constructively, because 
they are coherent. This is sometimes called first order interference. We 
will see in the following that the Bose-Einstein correlations of interest 
for us are due to the incoherence of the production process, and are a 
manifestation of the so called second order interference. 

• The statement that the scattering probability for identical particles is 
twice the corresponding probability for distinguishable particles is not 
possible to check experimentally, because non identical alpha particles 
are not available. The best one can do is to compare the experimental 
result for the identical alpha particles with the calculation for the non 
identical ones. In the present example, where the calculation is simple 
and non controversial, this is not much of a problem, but in multiple 
production processes a calculation from first principles is not possi- 
ble and the definition of the distribution for distinguishable particles, 
which should be modified by the Bose-Einstein correlations to yield the 
distribution which can be compared with experiment, is a difficulty 

• The final state can be represented by the density operator 
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p = -\U X L 2 + U 2 L X ){U X L % + U 2 L X \, (1) 

where £/jL fc is the state, where particle % is registered by the detector U 
and particle k by detector L. Expanding the left-had-side one obtains 
four terms, if, however, the density operator is to be used only for 
calculating averages of operators symmetric with respect to exchanges 
of the identical particles, which is sufficient for all practical applications, 
one can use the simpler form 

p=\U x L 2 )(U x L 2 \ + \U x L 2 )(U*L 1 \. (2) 
It is useful to rewrite this formula in the form 

p = Y,\UiL 2 )(U P1 L P2 \, (3) 
p 

where the summation is over all the permutations P of the indices 1 and 
2 and Pi denotes the index obtained from index i under permutation 
P. For our simple example this formula is ridiculously complicated, for 
more difficult cases, however, its analogues are very convenient. 

2 HBT contribution 

An interesting application of the Bose-Einstein interference to find the sizes 
of the emitting objects was discovered in the fifties by two astronomers, R. 
Hanbury Brown and R.Q. Twiss By studying the second order Bose- 
Einstein interference of photons, they were able to measure the radii of some 
stars. The idea may seem obvious today, but it was not so at the time it was 
put forward. In the seventies Hanbury Brown wrote (quoted after ||) 

Now to a surprising number of people this idea seemed not only 
heretical but patently absurd and they told us so in person, by 
letter, in publications, and by actually doing experiments which 
claimed to show that we were wrong. At the most basic level they 
asked how, if photons are emitted at random in a thermal source 
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can they appear in pairs at the two detectors. At a more so- 
phisticated level the enraged physicist would brandish some sacred 
text, usually by Heitler, and point out that . . . our analysis was 
invalidated by the uncertainty relation . . . 

Since the distances to many stars are known, their radii could be deter- 
mined, if the opening angles between the light rays coming from the stars 
could be measured. These angles, however, are in most cases too small for 
a direct measurement. Hanbury Brown and Twiss suggested the following 
procedure. Consider two light rays coming from two points on the surface 
of the star - ray a from point a and ray b from point b. The problem is to 
measure the angle between the two rays. Each of the rays falls on two 
photodetectors denoted 1 and 2. The distance between the photodetectors 
is d. Elementary trigonometry yields to first order in the relation 

e = (4) 

a sin a 

where A$ is the difference of distances between point % on the star and the 
two photodetectors, while a is the angle between the line connecting the 
two photodetectors and the direction of the two rays. Thus the problem of 
measuring the opening angle reduces to the problem of measuring A a — 
A b . Of course in practice, in order to find the radius of the star a suitable 
averaging over the possible emission points is necessary, this, however, is 
rather simple and we shall not discuss it any further. 

The current generated in the photodetector is proportional to the inten- 
sity of the incident light. Thus for counter 1 it is 

i lu = K x [E a sin(cj a t + fa) + E b s\n(u h t + fa)] 2 , (5) 

where fa, E iy oji denote respectively the phase at the star surface, the am- 
plitude and the frequency for ray i, and K\ is a proportionality coefficient 
dependent on the working of the photodetector 1. For simplicity, the polariza- 
tion effects have been ignored and the time necessary to reach photodetector 
1 has been put equal t for both point a and point b. In the apparatus the 
current i\ u is further filtered so that only frequencies from 1 Hz to 100 Hz 
survive. Thus finally, the current from the first photodetector is 

ii = K x E a E b cos[(w a - Lu b )t + (fa - fa)]. (6) 



This current is zero on the average and does not look particularly interesting. 
The analysis for the second photodetector is similar except that the time 
necessary to reach the detector for the ray from point % on the star is increased 
by Aj/c . One obtains 



i 2 = K 2 E a E b cos 



{u a -u; b )t + -(A a - A b ) + 
c 



(7) 



where uo « u a ~ 0J b . This is again a rather uninteresting current, but the 
average of the product of the filtered currents from the two photodetectors 



(iii 2 ) = K x K 2 E a E b cos 



OJ 



(A 6 - A 



(8) 



which is measurable, yields A a — A b and consequently the necessary opening 
angle 0. 

Note that the result is obtained in spite of the fact that the presence of the 
random phases 0j means that light from a is incoherent with respect to light 
from b. Because of these phases the product of two amplitudes, one for the 
ray a and one for the ray b averages to zero, The product of four amplitudes, 
two from a and two from b, however, can survive. For this reason one calls 
this effect second order interference or intensity interferometry. 



3 The GGLP contribution 

The first application of intensity interferometry in particle physics was made 
by the Goldhabers Lee and Pais ||. Their problem was somewhat different 
from that that in the HBT case. The interfering particles were like sign pion 
produced at two points within the interaction region of a hadron - hadron 
collision. The interference of interest was not between the measurements 
at two points in space, but between momentum measurements. Assuming 
that at the production space-time point x^ the pion wave function has phase 
0^ and that the momentum of the pion p^ is well-defined, one expects at 
the registration point x an amplitude proportional to exp[ipk(xk — x) + i<frk\- 
The probability of finding the two pions produced at points x\ and x 2 with 
momenta pi and p 2 , after proper symmetrization of the wave function, is 
proportional to 
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= 1 + cos[(pi -p 2 )(xi - x 2 )]. (9) 

We assume now that the production process is incoherent, so that the av- 
eraging over the times and positions xi, x 2 should be made at the level of 
probabilities and not of amplitudes. Then the distribution of the difference 
in momenta should be approximately given by the formula 

C(Pi -P2) = 1+ (cos[(pi -p 2 ){xi - x 2 ))}, (10) 

where the averaging is over X\ and x 2 . Qualitatively, the result does not 
depend much on the actual prescription being used for the averaging. For 
p 1 m p 2 the argument of the cosine is close to zero and consequently C ~ 2. 
For large momentum differences, the argument of the cosine is a rapidly 
oscillating function of X\ — x 2 , which is strongly suppressed by the averaging 
process, and C ~ 1. If the weight function used for the averaging contains 
just one parameter with the dimension of length, let us denote it R, the width 
of region in p\ — p 2 , where C is significantly bigger than one, must be of the 
order of R" 1 . There are many specific recipes how to perform the averaging. 
The results obtained for the correlations of momenta and for the sizes and 
shapes of the interaction regions are reasonable. For reviews see f| and |J. 
In spite of this success many difficulties remain. 

• Since the energy of a pion is determined by its momentum, one has 
data only on the three-dimensional distribution of the differences of 
spacial momenta. This is not enough to derive the four-dimensional 
distribution of the sources in space-time. Therefore, the results are 
strongly model dependent. 

• The averaging over the square of the wave function corresponds to the 
assumption that the density matrix of the final pions in coordinate 
representation is diagonal. This in turn implies that the momentum 
distribution should be flat, which contradicts experiment. A closely 
related question is, how the pion can be initially localized at the pro- 
duction point and then represented by a plane wave corresponding to 
well-defined momentum. 

• Information about the production region is, in practice, obtained only 
from pairs of pions with similar momenta. Consequently, what is being 



e i(piXi+p2X2) _|_ e i(pix 2 +P2Xl) 
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measured is not the whole interaction region, but the region, where the 
pions with similar momenta are produced. This, incidentally, explains 
the fact that the interaction regions usually come out roughly spherical, 
while one believes that the full interaction region is more string like. 

• There are many corrections, which probably should be applied, but it 
is controversial how. Here belong the corrections for coulomb repulsion 
between the charged like sign pions, the corrections for the final state 
interactions due to strong coupling, the corrections due to resonance 
production, the corrections due to partial coherence of the source etc. 



4 Density matrix approach 

In order to obtain a more general formulation for the GGLP problem it is 
convenient to use the formalism of density matrices. This has been described 
by a number of people, here we are the closest to the formulation used by 
Bialas and Krzywicki j7|. We introduce an auxiliary, unphysical process, 
where all the particles produced are distinguishable. We assume that for this 
process simple, intuitive ideas work. Then we correct for the Bose-Einstein 
correlations in order to obtain results comparable with experiment. This 
approach has its defects as discussed in the Introduction (see also [|], ||), 
but for lack of a better idea it is widely used. As our starting point for 
the distinguishable particles we use an independent production model (cf 
0) 0i E3 an d references contained there). In this model the multiplicity 



distribution for the particles is poissonian 

and for each multiplicity the density matrix is a product of single particle 
density matrices 

TV 

p£W) = npi 0) (^)- (12) 

i=i 

The momentum distribution is given, as usual, by the diagonal elements of 
the density matrix in the momentum representation 
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n 0N (q) = Pn (q,q). (13) 
It is convenient to normalize it to unity 

aP = JdqQ 0N (q) = l. (14) 

For identical particles the density matrix should be symmetrized as ex- 
plained in the Introduction 

pnM) = T,p*Mp)- (is) 

p 

The corresponding momentum distribution for a given multiplicity is 

^N = PN(q,q)- (16) 
This, however, is no more normalized to unity, because 

a N = J dqtt N {q) = 1 + ... (17) 

The first term in the last expression corresponds to the identity permutation, 
but there are (N — 1)! further terms. This yields the multiplicity distribution 

P N =tfP$ ) a N , (18) 

where M is an iV-independent normalizing factor, which ensures that J] Pn = 
1. 



5 Simple case: pure final state 

In order to present simply the qualitative features of the result, let us consider 
first the case, when for each multiplicity the final state is pure. This is a 
grossly oversimplified model, but we will find that it contains some features 
of the much more realistic approach presented in the following section. For 
the pure state model 

^ = \< ] )m (19) 
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We assume that each of the state vectors is symmetric with respect 

to exchanges of particles. Thus the effect of the summation over the per- 
mutations P is simply to multiply the operator p$ by N\. As a result the 
probability of producing exactly N particles is also multiplied by the factor 
N\. The Poisson distribution goes over into a geometrical distribution and 
after evaluating the normalization factor we get 

P N = {1- v)v N . (20) 

This formula makes sense only if v < 1, because otherwise the sum of the 
probabilities Pn diverges. For the average number of particles one finds 

37 = < 21 > 

with a singularity at v — 1. From the model presented in the following section 
it will be seen that this singularity corresponds to Einstein's condensation. 

In order to avoid the repeated summation of series it is convenient to 
introduce the generating functions. The generating function for the multi- 
plicity distribution is 

°° 1—7/ 

*(z) = £ Pnz n = —— (22) 

The logarithmic derivative of this function with respect to z at z — 1 yields 
the average multiplicity. The second derivative of the logarithm with respect 
to z at z = 1 is the dispersion and in general the p-th cumulant of the 
multiplicity distribution is given by 

Inclusive and exclusive momentum distributions, as well as all the corre- 
lation functions in momentum space, can be calculated by functional differ- 
entiation from the generating functional 



oo „ N i 

$M = £ MP^ / dqn N (q) Hu( qi ) = (24) 

For instance, the single particle distribution is 
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\ / M=l 

Thus symmetrization (Bose-Einstein statistics) introduces in this simple model 
only a change of normalization. 

6 Independent production 

Let us consider now the full independent production model. In order to find 
the modification of the multiplicity distribution due to Bose-Einsten statistics 
it is necessary to calculate the correction factors 

r N 

a N = / d( l lM 0) (ft>9«)- ( 26 ) 
p J 1=1 

Since each permutation can be decomposed into cycles, this integrals can be 
expressed in terms of the cycle integrals 

C k>1 = J d 3k q pf\ qi , q2 )pf\ q2 ,q 3 )...pf\ qk , qi ). (27) 
It is convenient to add the definition 

d = 1. (28) 

Similarly the integrals necessary to calculate the generating functional for 
the momentum distributions can be expressed in terms of the cycle integrals 

C k [u] = J d 3k q u( qi )p?\ qi , q2 )u( q2 )p?\ q2 , q3 ) . ..u(Qk)A%k,Qi) (29) 

After some combinatorics, very similar to that used when deriving the 
linked clusters expansion familiar from quantum field theory and many body 
theory, one finds the generating functional 



$[it] = exp 



- ^[u] - C k [l] 

" h 



(30) 



Substituting in this functional z for u one obtains the generating function 
for the multiplicity distribution. Without exhibiting the actual calculations 
we will now present some general results, obtained for this model. 
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• The single particle momentum distribution and all the momentum cor- 
relation functions can be expressed in terms of one function depending 
on two single particle momenta 



HquQi) = / d 3 q 2 ...d 3 q k pY ) (qi,q2)pi(q2,q3)---Pi (qk,q' 1 )- 

k=l J 

(31) 

For instance the momentum distribution is 

n(q) = L(q,q). (32) 
The two particle cumulant is 

K 2 (qi, q 2 ) = L(q u q 2 )L(q 2 , qi). (33) 
In general the p-th correlation function is 



K p (qi, q 2 ,..., q P ) = L(q u q 2 )L(q 2 , q 3 ) . . . L(q p , q ± ) + 

(permutations of the indices 2, . . . ,p). (34) 

• For typical density matrices the average square of the difference of 
momenta between two particles (q 2 ) decreases due to symmetrization. 

• For typical density matrices the average difference between the produc- 
tion points of pairs of particles decreases due to symmetrization. 

• For typical density matrices the size of the interaction region as evalu- 
ated from the width of the two-particle momentum correlation function 
decreases due to symmetrization. 

We will discuss these predictions in a further section, where we will rederive 
them in a more intuitive way. The references to "typical density matrices" 
mean that the statement is true for most density matrices, but not for all. We 
have not been able to find a condition defining the relevant class of density 
matrices. 

Probably the most interesting implication is the possibility of Einstein's 
condensation, but this will be discussed in the following section. 
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7 Einstein's condensation 



Using matrix notation one can rewrite the definition of the function L given 
in the previous section in the form 

i(?,9') = E^<9l(M 0) ) W)- (35) 
k=i 

Expanding the single particle density operator in terms of its eigenvectors 
and eigenvalues, we find 

= £ \ n )X n (n\ (36) 

n 

and for its k-ih power 

(p! 0) ) fe = Y,MH 07) 

n 

Thus in the momentum representation 

00 

L(g,gO=E<^>W>E^- (38) 

n k=l 

This expression makes sense only if for all n there is \ n v < 1. Denoting the 
largest eigenvalue of the density operator by Ao, we expect problems when 
u\ ->• 1. 

Performing the summations of the geometric series, we can rewrite the 
expression for L(q, q') in the form 

L( q , g ') = Y, Uqm( ( )UXn - (39) 
For v\ — > 1 it is convenient to use the equivalent formula 

LM)= Mimm +i(q , q% (40) 

1 — VAq 

where L remains bounded in the limit. Putting q — q' and integrating over 
q we get the corresponding formula for the average multiplicity 

AT = - — _|_ bounded term. (41) 

1 ~ v A n 
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From these formulae it is clear that when z/Ao tends to one, Einstein's con- 
densation occurs. Increasing v corresponds to the increasing of the number 
of particles in the system. In all the states with indices n^O there is place 
only for a limited number of particles, while all the surplus, which can be 
arbitrarily large, gets located in the state |0). In this sense, when the number 
of particles becomes very large, we recover the model with the pure state dis- 
cussed previously. A very interesting question is, whether experimentally it is 
possible to create condition, where the Einstein condensate would dominate. 

Let us conclude this section with two remarks. For the Gaussian single 
particle density matrix 



(0) / t\ 

Pi {q,q) = /— r^ exp 

V ITILX 



"IT 9- 



(42) 



2A 2 2 

where q + = (q + q')/2 and q~ = q — q', the eigenvalues and eigenfunction are 
known ||. Thus all the calculations can be easily performed. In fact they 
have been performed by various methods JTT]], flI2| , [T3|], |TJJ. 

The theory can be reformulated in the second quantization formalism. 
Then the function L(q,q') appears as the Green function {a) q a q i) and the 
possibility of expressing all the correlation functions in terms of L(q, q') is 
the Wick theorem with L as the only non zero contraction. 



8 Statistical physics interpretation 

Many results from the previous sections can be simply reinterpreted and 
rederived using standard statistical physics. Consider the single particle un- 
symmetrized density operator 

pS = EI">M«I ( 43 ) 

n 

with the condition J2n — 1- This can be reinterpreted as the density 
operator corresponding to the canonical ensemble, if we put 

An = |e"^, (44) 

where, as usual, e n is the energy of state |n), (3 is the inverse temperature in 
energy units and Z = J2 n exp[— j3e n ] is the canonical partition function. The 
corresponding (single particle) Hamiltonian is 
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if = X»e n (ra|. (45) 

n 

This Hamiltonian, when written in the coordinate representation, may look 
quite unusual, but some cases are simple. For instance, the Gaussian density 
matrix corresponds to the Hamiltonian of a harmonic oscillator. 

For indistinguishable particles a single particle is not a convenient subsys- 
tem and, as suggested by Pauli long ago, it is better to choose as subsystem 
the open system consisting of all the particles in state \n). The state of this 
subsystem is defined by the number of particles (N) in it. The probability 
of state N of the subsystem is 

P n (N) = —v N e- pN£n . (46) 
The grand partition function 

Z » = 1 - -Be ( 4? ) 

1 — ve pe 

is chosen so that J2n=o Pn(N) = 1, the parameter 17 is known in statistical 
physics as the fugacity and is connected to the chemical potential /x by the 
formula 

V = e^. (48) 

In order to reproduce the formulae from the previous sections, one puts 
Zv = v. The grand partition function can be used to find the moments 
of the multiplicity distribution very much like the multiplicity generating 
function. For instance, for the average occupation of state n we find 

_l dLogZ n 1 v\ n 

™ dfx e ^»-f)-i 1-i/V 1 yj 

The probability of no particles in state n is 

The probability of no particle in the whole system is 
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' (0) =n<*jTT- (51) 

Let us consider two limiting cases. When all the occupation numbers are 
very small, the product equals approximately exp[— (N)] and for large mul- 
tiplicities it is very small. Very large fluctuations of the multiplicity are very 
unlikely to occur. When most particles are in the state n — 0, the product is 
approximately ((N) + which is much bigger than in the previous case. 
Thus, when there is much Einstein condensate, large multiplcity fluctuations 
become much more probable. Cosmic ray physicists have been reporting f[5 



observations of centauro and anticentauro events. This are high multiplcity 
events, where respectively either the neutral pions or the charged pions are 
missing. One could speculate that this phenomena are related to Einstein's 
condensation. 

Statistical physics gives also a simple interpretation for the function L(q, q') 
One finds 

L(q,q') = Ei(^(?l^)7T' (52) 

This is the canonical density matrix with the Maxwell-Boltzmann weights re- 
placed by the Bose-Einstein weights. The fact that the Bose-Einstein weights 
fall with increasing energy e n faster than the Maxwell-Botzmann weights ex- 
plains qualitatively most of the observations reported previously. For most 
Hamiltonians the wave function spreads in ordinary space and in momentum 
space, when energy is increased. Since the Bose-Einstein weights enhance the 
low energies, they reduce the average momenta and radii. Also the reduction 
of the effective radius of the interaction region, as determined from the width 
of the correlation function in momentum space, can be easily understood. If 
in the previous formula all the terms had equal weights, we would obtain 
L(q, q') = 5 3 (q — q'). The stronger the cut on the sum, the broader the peak 
in q — q' becomes. Since the Bose-Einstein weights are more peaked at low 
energies than the Maxwell-Boltzmann ones, they correspond to a broader 
peak in the correlation function. Since the width of this peak is inversely 
proportional to the radius of the production region, symmetrization reduces 
the radius of this region. All these qualitative arguments are usually true. 
It is, however, easy to show examples of hamiltonians, where e.g. with in- 
creasing energy the wave function shrinks either in ordinary space, or in 
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momentum space. Additional assumptions necessary to convert these qual- 
itative arguments into rigorous theorems are, therefore, necessary, but not 
yet known. 
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